• Multi-site imaging studies are becoming increasingly common.
• Combining imaging data across sites introduces non-biological sources of variation that arise from the use of different scanner hardware and acquisition protocols.
  • E.g., field strength, manufacturer, subject positioning
• Scanner effects or site effects are similar to batch effects in the genomics literature
  • Known to affect measurement of regional volumes, cortical thickness, voxel-based morphometry, …
  • More generally, structural, functional, diffusion tensor, and other types of images and features extracted from them may exhibit scanner effects

• Need to eliminate or account for scanner effects in downstream statistical analyses
  • Most critical if sites or scanners are imbalanced with respect to other variables such as age, sex, race, clinical status
  • Simply including scanner as a confounding variable may not work well (Rao et al. 2017)
• Several methods for estimating and removing unwanted sources of variation due to site/scanner have been adapted to neuroimaging data.
• In this tutorial we will use ComBat (Johnson, Li, and Rabinovic 2007; Fortin et al. 2018) to harmonize cortical thicknesses from the ADNI data.

• Template: publicly available vs. "cooking" one based on your data
• DiReCT: Diffeomorphic Registration based CT measurement (Das et al. 2009)

• With tissue priors and N4 corrected T1 from atropos, can use cort_thickness in extrantsr package

• The Alzheimer's Disease Neuroimaging Initiative (ADNI) is a multi-million dollar study funded by public and private sources.
  • National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering, the Food and Drug Administration, private pharmaceutical companies, and non-profit organizations.
• The goals of the ADNI are to better understand progression from normal aging to mild cognitive impairment (MCI) and early Alzheimer's disease (AD) and determine effective biomarkers for disease diagnosis, monitoring, and treatment development.
• We estimated cortical thicknesses from a subset of initial subject visits (N=156)
  • Mix of male/female and healthy controls, MCI, and AD diagnoses at baseline
• Our subset consists of images aquired from scanners at 14 different sites
  • Mix of field strength, manufacturer, and model

• Multi-atlas label fusion applied using 20 OASIS template T1s manually labeled with the Desikan-Killiany-Tourville (DKT) cortical labeling protocol (www.mindboggle.info).

   label                               roi
1  X1002  “left caudal anterior cingulate”
2  X1003      “left caudal middle frontal”
3  X1005                     “left cuneus”
4  X1008          “left inferior parietal”
5  X1017                “left paracentral”
6  X1022                “left postcentral”
7  X1023        “left posterior cingulate”
8  X1024                 “left precentral”
9  X1025                  “left precuneus”
10 X1027     “left rostral middle frontal”
11 X1028           “left superior frontal”
12 X1029          “left superior parietal”
13 X1030          “left superior temporal”
14 X1031              “left supramarginal”
15 X1035                     “left insula”

• Scanner is often confounded with biological covariates of interest.

• ComBat (Johnson, Li, and Rabinovic 2007) has been shown to remove unwanted variation due to scanner while preserving biological associations (Fortin et al. 2018).
• ComBat assumes expected values of the imaging feature measurements can be modeled as a linear combination of the biological variables and the site effects with an error term that includes a multiplicative scanner-specific scaling factor: \[y_{i,j,v} = \alpha_v + X_{i,j}^{T}\beta_v + Z_{i,j}^{T}\theta_v + \delta_{i,j}\epsilon_{i,j,v}\]
  • \(y_{i,j,v}\) is average cortical thickness in the \(v^{th}\) ROI from subject \(i\), scanner \(j\)
  • \(X_{i,j}^{T}\) is a vector of fixed covariates for subject \(i\) scanned on scanner \(j\)
  • \(\alpha_v\), \(\beta_v\) are the intercept and slope vector of covariates for the \(v^{th}\) imaging measurement
  • \(Z_{i,j}^{T}\) is a vector that indicates which scanner subject \(i\) was imaged on
  • \(\theta_v\) is the vector of additive scanner effects
  • \(E(\epsilon_{i,j,v}) = 0\)

• Substituting \(\gamma_{i,v}\) for \(Z_{i,j}^{T}\theta_v\) we have: \[y_{i,j,v} = \alpha_v + X_{i,j}^{T}\beta_v + \gamma_{i,v} + \delta_{i,j}\epsilon_{i,j,v}\]
• ComBat uses empirical Bayes estimation to improve quality of scanner-specific parameter estimates when sample sizes are small
  • Let \(\gamma_{i,v}^{*}\) and \(\delta_{i,v}^{*}\) denote the EB estimates of \(\gamma_{i,v}\) and \(\delta_{i,v}\), respectively.
• Then, the ComBat-harmonized cortical thicknesses are defined as \[y_{i,j,v}^{ComBat} = \frac{y_{i,j,v} - \hat{\alpha}_v - X_{i,j}^{T}\hat{\beta}_v - \gamma_{i,v}^{*}}{\delta_{i,v}^{*}} + \hat{\alpha}_v + X_{i,j}^{T}\hat{\beta}_v\]

• Need utils.R and combat.R from https://github.com/Jfortin1/ComBatHarmonization

source('utils.R')
source('combat.R')

• The model matrix should contain covariates of biological interest.

modelData = read.csv('modelData.csv')
head(modelData)

     subject     age sex     dx            scanner
1 002_S_0413 76.4329   F Normal 3_PHILIPS_Intera_2
2 002_S_0559 79.4658   M Normal 3_PHILIPS_Intera_2
3 002_S_0729 65.2986   F    MCI 3_PHILIPS_Intera_2
4 002_S_0816 70.9534   M     AD 3_PHILIPS_Intera_2
5 002_S_0954 69.5041   F    MCI 3_PHILIPS_Intera_2
6 002_S_1018 70.8055   F     AD 3_PHILIPS_Intera_2

• We will include age, sex, and diagnosis.

mod = model.matrix(~age+factor(sex)+factor(dx), data=modelData)

• Let's read in the cortical thickness data.

ctData = read.csv('imageData.csv')
head(ctData)[,1:10]

     subject    X1002    X1003     X1005    X1008    X1017    X1022
1 002_S_0413 2.349515 1.957340 1.6771022 2.929418 1.893073 1.733991
2 002_S_0559 2.814481 1.768518 0.9173002 2.297948 1.612640 2.071636
3 002_S_0729 2.788202 2.108902 1.6343228 3.154604 2.219242 1.984391
4 002_S_0816 2.753477 2.168249 1.7955870 2.880065 1.973897 2.042263
5 002_S_0954 2.523273 1.664635 1.0189855 2.077749 1.236037 1.468252
6 002_S_1018 2.596688 2.216399 1.9206076 2.424231 1.744081 1.800574
     X1023    X1024    X1025
1 2.365209 1.968747 2.474627
2 2.606214 1.838712 2.345573
3 2.578613 1.782150 2.055732
4 2.188748 1.990839 2.992091
5 1.329828 1.468782 1.623318
6 2.165901 1.947617 2.611285

• The imaging data needs to be in a separate matrix where rows are features and columns are subjects.
• Let's remove the subject column and transpose the ctData object.

img = t(ctData[,-1])
head(img)[,1:10]

          [,1]      [,2]     [,3]     [,4]     [,5]     [,6]     [,7]
X1002 2.349515 2.8144805 2.788202 2.753477 2.523273 2.596688 1.833424
X1003 1.957340 1.7685180 2.108902 2.168249 1.664635 2.216399 1.696198
X1005 1.677102 0.9173002 1.634323 1.795587 1.018986 1.920608 1.228951
X1008 2.929418 2.2979484 3.154604 2.880065 2.077749 2.424231 2.712218
X1017 1.893073 1.6126404 2.219242 1.973897 1.236037 1.744081 1.760396
X1022 1.733991 2.0716363 1.984391 2.042263 1.468252 1.800574 1.712189
          [,8]     [,9]     [,10]
X1002 2.960386 3.191757 1.6549346
X1003 2.220978 2.294089 1.1845660
X1005 1.638375 1.760091 0.6901224
X1008 2.727885 2.280574 1.2541459
X1017 2.105533 1.890207 1.0815922
X1022 2.095173 1.480126 0.8755437

• Given a model matrix, scanner information, and image data matrix, ComBat just requires a single line of code:

harmonized = combat(dat=img, batch=modelData$scanner, mod=mod)

[combat] Performing ComBat with empirical Bayes
[combat] Found 14 batches
[combat] Adjusting for 4 covariate(s) or covariate level(s)
[combat] Standardizing Data across features
[combat] Fitting L/S model and finding priors
[combat] Finding parametric adjustments
[combat] Adjusting the Data

• 4 covariates: age, sex (2-level factor), and diagnosis (3-level factor).
• The combat() function returns a list.

• The combat() function returns a list.
• Harmonized data are returned as a matrix with the same dimensions as the image data input (rows are features, columns are subjects).

head(harmonized$dat.combat)[,1:10]

          [,1]      [,2]     [,3]     [,4]      [,5]     [,6]      [,7]
X1002 2.075969 2.5403388 2.514211 2.479313 2.2494910 2.322785 1.5601570
X1003 1.604147 1.3928650 1.774367 1.829854 1.2947712 1.890831 1.3247258
X1005 1.251037 0.6043039 1.215196 1.363586 0.6807333 1.459647 0.8726338
X1008 2.242658 1.6391016 2.457074 2.191413 1.4269869 1.755756 2.0304396
X1017 1.414838 1.1241720 1.749303 1.497199 0.7381837 1.263199 1.2756135
X1022 1.172603 1.5470296 1.459725 1.521212 0.8741849 1.260857 1.1429860
          [,8]     [,9]     [,10]
X1002 2.686305 3.134451 1.7216344
X1003 1.886106 2.388280 1.3679891
X1005 1.219925 1.922309 0.8532241
X1008 2.053663 2.445597 1.4071359
X1017 1.632919 2.067534 1.2650911
X1022 1.579223 1.739706 1.1057018

• Compare prior distributions to EB estimated parameter distributions
  • Location: harmonized$gamma.hat versus harmonized$gamma.star
  • Scale: harmonized$delta.hat versus harmonized$delta.star

• Before ComBat

• After ComBat

• Before ComBat

• After ComBat

• Before ComBat

• After ComBat

modelData$sex = factor(modelData$sex)
modelData$dx = factor(modelData$dx)
modelData$scanner = factor(modelData$scanner)

preComBat = left_join(modelData, ctData, by='subject')
preRInsula = lm(X2035 ~ age + sex + dx + scanner, data=preComBat)
summary(aov(preRInsula))

             Df Sum Sq Mean Sq F value   Pr(>F)    
age           1  4.222   4.222  24.157 2.48e-06 ***
sex           1  0.297   0.297   1.701  0.19427    
dx            2  2.007   1.004   5.743  0.00402 ** 
scanner      13 12.456   0.958   5.482 6.07e-08 ***
Residuals   138 24.119   0.175                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

harmonizedData = as.data.frame(t(harmonized$dat.combat))
harmonizedData$subject = ctData$subject

postComBat = left_join(modelData, harmonizedData, by='subject')
postRInsula = lm(X2035 ~ age + sex + dx + scanner, data=postComBat)
summary(aov(postRInsula))

             Df Sum Sq Mean Sq F value  Pr(>F)   
age           1  1.625  1.6250  10.377 0.00159 **
sex           1  0.550  0.5497   3.511 0.06309 . 
dx            2  1.799  0.8994   5.744 0.00402 **
scanner      13  0.529  0.0407   0.260 0.99566   
Residuals   138 21.610  0.1566                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Data harmonization is an important step in any image analysis that combines data from different sites and/or scanners.
• ComBat has been successfully applied to DTI and cortical thickness data to remove scanner effects while preserving biological associations of interest (Fortin et al. 2017; Fortin et al. 2018).
  • ComBat code is located at https://github.com/Jfortin1/ComBatHarmonization
  • Available in R and MATLAB
• To obtain valid statistical inference in downstream analyses, uncertainty in ComBat harmonization should be considered

http://johnmuschelli.com/imaging_in_r